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' Abstract 

' The Westervelt equation is a model for the propagation of finite 

d . amplitude ultrasound. The method of discrete exterior calculus can be 

used to solve this equation numerically. A significant advantage of this 
method is that it can be used to find numerical solutions in the discrete 
space manifold and the time, and therefore is a generation of finite 
difference time domain method. This algorithm has been implemented 
(N| 1 in C++. 
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1 Introduction 

•l-H 

j_j ■ The nonlinear effects are important in medical ultrasound and also plays a 



role in diagnostic program. Some of diagnostic ultrasound instruments have 
implemented harmonic imaging into their devices by receiving the harmon- 
ics in the reflected ultrasound caused by nonlinear distortion of the signal 
propagating through the biological tissue. Meanwhile, we should pay atten- 
tion to observing hygienic limits of ultrasound exposition to avoid heating 
the unwanted tissue. The ability to predict the effects of nonlinear ultra- 
sound propagation therefore becomes important. Numerical simulations are 
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currently the best means of making predictions of nonlinear ultrasound prop- 
agation. Two of the numerical techniques suitable for providing solutions 
to the propagation problem are the finite difference time domain (FDTD) 
method [IH3] and the discrete exterior calculus(DEC) [3HT2]. 

The Westervelt equation is a model for the propagation of finite ampli- 
tude ultrasound, deriving from the equation of fluid motion by keeping up to 
quadratic order terms |13}ll4j. In this paper, an conditional stable scheme 
for this equation using the techniques in DEC is proposed. A significant 
advantage of this scheme is that it can be used to find numerical solutions 
on the discrete space manifold and the time, therefore is a generation of 
FDTD. In order to refine the space mesh to improve the accuracy of the 
solution given by the this scheme, the amount of work involved increases 
rapidly, since the length of time step should also be reduced. Hence, two 
unconditional stable schemes are also proposed, namely the implicit and 
semi-implicit DEC schemes for Westervelt equation on manifold. 

2 DEC schemes for Westervelt equation 
Discrete Laplace operator 

The 2D or 3D space manifold can be approximated by triangles or tetra- 
hedrons, and the time by line segments. Suppose each simplex contains its 
circumcenter. The circumcentric dual cell D{a$) of simplex <ro is 

D(a ) := [J Int(c(cr )c(cJi) • • • c(oy)), 

<xoGai6---£oy 

where crj is all the simplices which contains (To,..., 0i-i, and c(<7j) is the 
circumcenter of <jj. A discrete differential A:-form, k G Z, is the evaluation of 
the differential fc-form on all fc-simplices. Dual forms, i.e., forms evaluated 
on the dual cell. In DEC, the exterior derivative d is approximated as the 
transpose of the incidence matrix of fe-cells on k + 1-cells, the approximated 
Hodge Star * scales the cells by the volumes of the corresponding dual and 
primal cells, and the Laplace operator is approximated as 

A w *~ l d T * +d T * d. 

Take Fig.l as an example for a part of 2D mesh, in which 0,..., C are vertices, 
1, 2, 3 are the circumcenters of triangles, a, b, c are the circumcenters of 
edges. Denote kj as the length of line segment and Ayjy as the area of 
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quadrangle (i,j,k,l). 




Define 
and 



Figure. 1. A part of 2D mesh 

ll2 '■= hb + hb, ^23 : = he + he, hi '■= ha + ha, 



Pl23 '■— Aoi a b + A()2bc + A)3ac- 

In Fig.l, the discrete Laplace operator acting on p at vertex 

1 / l\3 l\2 / « ^23 

A P0 ~ "75 -j—{PA - PO) + 1 [PB ~ PO) + -j—{PC- 

"123 \'A0 'BO 'CO 





(I.pj 







Figure. 2. A part of 2D rectangular mesh 
In Fig. 2, scheme (1) reduces to 

A Phj+l +Pi,j-l +Pi+l,j +Pi-l,j ~ tyij 

Pi ' j ~ (As) 2 
where As is the uniform space step length. 

Discretization for Westervelt equation 

The Weservelt equation can be written in the following form: 

1 d 2 p 5 d 3 p f3 d 2 p 2 

p ~c 2 W 2+ 7 W + p 4 dt 2 = ' 



where p is the acoustic pressure, po is the are the ambient density, and cq is 
the ultrasound speed, 5 is the diffusivity of ultrasound, j3 is the coefficient 
of nonlinearity. 

For some situations, a source having azimuthal symmetry about its axis 
is considered. In this case, 2D triangular discrete manifold as the space is 
only need to be considered. And the Eq.(2) on 3D space manifold and the 
time can also be solved numerically, using a similar approach. 

The temporal partial derivatives with discrete differences, which can be 
obtained from Taylor series expansions about each node of the computa- 
tional mesh. Temporal derivatives present in Eq.(2) can be calculated as 
follows: 

^(A^" +1 -V +*>""') 

x3 n -1 

~ ;{p n - + 3 P n ~ 2 - P n - 3 ) 



Of 3 _ (At) 3 

.((p2 ) n_ 2(p 2 ) n-l + (p 2 



9 2 (p 2 ) n _ 1 /,2\n _ (^\n-l _i_ (J2\n-2\ 



dt 2 (At) 2 

where At is the time step length, n is the time coordinate. Let 



-Po'') - ^^((Po)"- 1 - ^pIY' 2 + (PIT- 3 ) 



4(At)3^0 u ^0 

/,,,-'>» I _ 

The explicit DEC scheme for Eq.(2) is 
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= Lpg" 1 . (3) 



If refining the space mesh, the amount of work involved increases rapidly, 
since the length of time step should satisfy the stable condition. Now, two 
unconditional stable schemes are proposed, namely the implicit DEC scheme 
(4) and semi- imp licit DEC scheme (5). 

^- (r-iPl-PZ) + i^-(pb-po) + r^(P n c-p n )) = Lpr 1 (4) 

-n.23 VAO I BO 'CO / 

The higher order accuracy scheme for temporal derivatives can also be 
used in schemes(3-5). 
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3 Stability, convergence and accuracy 



The first two terms in (2) describe linear wave propagation at the small- 
signal sound speed. The third term describes the loss due to he viscosity 
and thermal conduction of the media. The fourth term describes the non- 
linear distortion of the traveling wave due to amplitude effects. Since the 
coefficients of third and fourth terms in Eq.(2) are small compared with the 
coefficient of second term, the effect of ^§ and ^Jr on the stability of the 
scheme (3) can be ignored. Therefore, the stable condition for the wave 
equation 

1 d 2 p 

is need to be satisfied. The wave equation can be discretized as 

RightUr- 1 = ^±_( p n-2p n - l +p" - 2 ) (6) 
The stable condition condition for scheme (6) [12] is 

C At < minvoevertices 



^ ^13 , _^12_ ^23_ 

\ P123 vao Ibo Icq 



(7) 



In Fig.2, the inequality (7) reduces to 

\/2 . 
c At < —As, 
~ 2 

which is the stable condition for scheme (6) on square grid. The analysis 
of unconditional stability for schemes (4) and (5) can also be proved as the 
wave equations. 

By the definition of truncation error, the exact solutions of Westervelt 
equation satisfy the same relation as schemes (3-5) except for an additional 
term O(At) 2 . This expresses the consistency, and so the convergence for 
schemes (3-5) by Lax equivalence theorem (consistency + stability = con- 
vergence). The derivative in Eq.(2) is approximated by first order difference 
in schemes (3-5). Equivalently, p is approximated by linear interpolation 
function. This is to say that schemes (3-5) have first order spacial and tem- 
poral accuracy, according to the definition about accuracy of finite volume 
method. 
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4 Algorithm Implementation 



The DEC schemes (3-5) for Westervelt equation was implemented in C++. 
The implementation consists of the following steps: 

1. Set the simulation parameters. These are the dimensions of the com- 
putational mesh and the size of the time step, etc.; 

2. Set the propagating media parameters. 

3. Initialize the mesh indexes. 

4. Assign current transmitted signal. 

5. Compute the value of all spatial nodes and temporarily store the result 
in the circular buffer for further computation. 

6. Visualize the currently computed grid of spatial nodes. 

7. Repeat the whole process from the step 4, until reach the desired total 
number of iterations. 

In the common practice, not every simulation step needs to be visualized, 
especially when the time step size is too small. In all the following examples, 
scheme (3) is used, and the parameters are 5 = 0.01 j3 = 1 p = 10000. The 
Gaussian envelope ^= exp(— y) cos(t) is used as a source signal. 

The example in Fig. 3 shows the nonlinear propagation effects and the 
diffusivity in the propagation of the Gaussian envelope. 



mmm 



Figure. 3. The propagation of Gaussian envelope on a rabbit 

The example in Fig. 4 shows the propagation of Gaussian envelope at 
a boundary of two kinds of media with different wave speeds 340m/s and 
3400m/s. 
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Figure. 4. The propagation of Gaussian envelope on a sphere with two 

kinds of media 
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